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ABSTRACT 

In order to study the process of cooling in dark-matter (DM) halos and assess how well simple 
models can represent it, we run a set of radiative SPH hydrodynamical simulations of isolated 
O ' halos, with gas sitting initially in hydrostatic equilibrium within Navarro-Frenk- White (NFW) 

potential wells. Simulations include radiative cooling and a scheme to convert high density 
t/^ , cold gas particles into collisionless stars, neglecting any astrophysical source of energy feed- 

^ I ' back. After having assessed the numerical stability of the simulations, we compare the result- 

ing evolution of the cooled mass with the predictions of the classical cooling model of White 
& Frenk and of the cooling model proposed in the MORGANA code of galaxy formation. We 
^ ' find that the classical model predicts fractions of cooled mass which, after about two central 

CO , cooling times, are about one order of magnitude smaller than those found in simulations. Al- 

though this difference decreases with time, after 8 central cooling times, when simulations are 
stopped, the difference still amounts to a factor of 2-3. We ascribe this difference to the lack 
of validity of the assumption that a mass shell takes one cooling time, as computed on the 
initial conditions, to cool to very low temperature. Indeed, we find from simulations that cool- 
ing SPH particles take most time in traveling, at roughly constant temperature and increasing 
density, from their initial position to a central cooling region, where they quickly cool down 
' to ~ 10"^ K. We show that in this case the total cooling time is shorter than that computed on 

the initial conditions, as a consequence of the stronger radiative losses associated to the higher 
, . ■ density experienced by these particles. As a consequence the mass cooling flow is stronger 

K> \ than that predicted by the classical model. 

The MORGANA model, which computes the cooling rate as an integral over the contribu- 
tion of cooling shells and does not make assumptions on the time needed by shells to reach 
very low temperature, better agrees with the cooled mass fraction found in the simulations, 
especially at early times, when the density profile of the cooling gas is shallow. With the 
addition of the simple assumption that the increase of the radius of the cooling region is coun- 
teracted by a shrinking at the sound speed, the MORGANA model is also able to reproduce 
for all simulations the evolution of the cooled mass fraction to within 20-50 per cent, thereby 
providing a substantial improvement with respect to the classical model. Finally, we provide 
a very simple fitting function which accurately reproduces the cooling flow for the first ~ 10 
central cooling times. 
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1 INTRODUCTION 

Understanding galaxy formation is one of the most important chal- 
lenges of modem cosmology. The rather stringent constraints on 
cosmological para meters now placed b y a number of independent 
observations (e.g., ISpringel et alJlOOd for a recent review) allows 
us to precisely set the initial conditions from which the formation 
of cosmic structures has started. As a consequence, understand- 



ing the complex astrophysical processes, related to the evolution 
of the baryonic component, represents now the missing link toward 
a successful description of galaxy formation and evolution. So far, 
two alternative approaches have been pursued to make quantita- 
tive predictions on the observational properties of the galaxy pop- 
ulation and their evolution in the cosmological context. The first 
one is based on the so-called semi-analytical models (SAMs, here- 
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after; e.g., Kauffmann et al 


19931: Somerville&Primackll 19991; 


ICole et al.l2000: Menci et al 


2005l;lMonaco et alj2007). In this an- 



proach, the background cosmological model predicts the hierarchi- 
cal build-up of the Dark Matter (DM) halos, where gas flows in, 
cools and gives rise to the formation of galaxies, while the complex 
interplay between gas cooling, star formation, chemical enrichment 
and release of energy from supernova (SN) explosions and Ac- 
tive Galactic Nuclei (AGN) is modeled through a set of simplified 
or phenomenological models, which are specified by a number of 
free parameters. A posteriori, the values of the relevant parameters 
should then be constrained by comparing SAM predictions to ob- 
servational data. The rather low computational cost of this approach 
makes it a quite flexible tool to explore the model parameter space. 

The second approach is based on resorting to full hydrody- 
namical simulations, which include the processes of gas cooling 
and suitable sub-resolution recipes for star formation and feedback. 
The obvious advantage of this method, with respect to SAMs, is 
that galaxy formation can be described by following in detail the 
gas-dynamical processes which determine the evolution of the cos- 
mic baryons during the shaping of the large-scale cosmic struc- 
tures. However, its limitation lies in the high computational cost, 
which makes it difficult to explore in detail the parameter space 
describing the process of galaxy formation and evolution. For this 
reason, following galaxy formation with full hydrodynamical sim- 
ulations in a cosmological environment of several tens of Mpc is 
a challenging t ask for simulations of the present generation (e.g.. 



Nagamine et al..2004; Na gai & Kravtsov 20 05; Romeo et al. 2005. ; 



Saro et al.ll200i 



This discussion shows that SAMs and hydrodynamical sim- 
ulations provide complementary approaches to the cosmological 
study of galaxy formation. The ability of hydrodynamical simu- 
lations of accurately following gas dynamics calls for the need of 
a close comparison between these two approaches, in order to test 
the basic assumptions of the SAMs. The best regime to perform 
this comparison is when one excludes the effect of all those pro- 
cesses, like feedback in energy and metals, whose different mod- 
eling in SAMs and simulation codes would make the comparison 
scarcely telling. Since gas cooling is the most basic ingre dient in 
any model of galaxy formation (e.g., IWhite & FrenkllT99lh . an in- 
teresting comparison would be performed when co oling is the only 
process turned on. In this spirit lBenson et aL lllooil), using a hydro- 
dynamical simulation of a cosmological box and a stripped-down 
version of SAM, compared the statistical properties of "galaxies" 
found in the two cases. They discovered that SPH simulation and 
SAM give similar results for the thermodynamical evolution of gas 
and that there is a very good agreement in terms of final fractions 
of hot, cold and uncollapsed gas. 

Similar conclusions were reached bv lHellv et all ( l2003h and 
ICattaneo et alj 1 2007 ). They improved the comparison performed 



by Benson et al. 1 200 ll) by giving to the down-stripped SAM the 



same halo merger histories extracted from the cosmological simu- 
lation. In this way they were able to compare cooling in DM halos 
not only statistically but on an object-by-object basis. The result 
was agai n that the two method s provide comparable "galaxy" pop- 
ulations. l^sUda^etal] (|20o2) performed a similar comparison for 
a simulation of a single galaxy cluster, obtaining similarly good 
results. 

While the general agreement between the two methods is en- 
couraging, still all the above analyses generally concentrated on 
comparing the statistical properties of the galaxy populations. Fur- 
thermore, if one wants to test the reliability of the cooling model 
implemented in the SAMs, the cleanest approach would be that of 



turning off the complications associated to the hierarchical merging 
of halos, thereby allowing gas to cool in isolated halos. 

The purpose of this paper is to present a detailed compari- 
son between the predictions of cooling models, as implemented in 
SAMs, and results of hydrodynamical simulations in which gas is 
allowed to cool inside isolated halos. Our controlled numerical ex- 
periments will be run for halos h aving the density profile (for DM 
particles) o f'Na varroetai] i ll 9971) (NFW hereafter), with a range of 
masses, concentration parameters and average densities (related to 
the halos' redshift). As a baseline model for gas c ooling, we con- 
sider the classical one, as originally proposed bv IWhite & FrenkI 
(1991). In this model, the cooling radius is defined as the radius 
at which the cooling time equals the time elapsed since radiative 
cooling is turned on. As a result, the growth rate of the cooled gas 
mass is simply related to the growth rate of th e cooling radius. This 
model was claimed bv lWhite&Frenkl ( ll99ll) to be close the exact 
self- similar solutions of cooling flows presented by Bertschinge^ 
( ll989D . Simulation results will also be compared to another model 
of gas cooling, which has been recently proposed bv Monaco et all 
( I2OO7I) in the context of the MORGANA SAM, and is based on a 
"dynamical" definition of cooling radius. As a result of our anal- 
ysis, we will show that the gas cooling rate in the simulations is 
initially faster than predicted by the classical cooling model. When 
the simulations are stopped, after about 8 central cooling times, this 
initial transient causes the classical cooling model to underestimate 
the cooled mass by an amount which can be as large as a factor 
of three, depending on the halo concentration, density and mass. A 
much better agreement with simulations is achieved by the alterna- 
tive MORGANA model. 

The plan of the paper is as follows. In sectio n 2 we describe 
first the "classical" analytic model for cooling (•white & FrenkI 
|l99 ll) . and the alternative MORGANA cooling model ( Monaco et al. 
2007). In section 3 we present numerical simulations, performed 
with the GADGET- 2 code and in section 4 we discuss the results 
obtained by comparing simulations to analytical cooling models. 
We discuss our results in section 5 and draw our final conclusions 
in section 6. A more technical discussion on the differences be- 
tween the classical and the MORGANA cooling models is provided 
in Appendix A, while Appendix B gives a very simple fitting for- 
mula for the cooling flows. 



2 ANALYTIC MODELS FOR GAS COOLING 

2.1 Hydrostatic equilibrium in an NFW halo 

All our tests start from a spherical dark matter halo with an NFW 
density profile, 



p{r) = Pcrit- 



(1) 



{r-/r,){l + r/rsY ' 

where is a scale radius, 5c is a characteristic (dimensionless) 
density, and pcrit is the critical cosmic density. The gas is assumed 
to be in hydrosta tic equilibrium. The equilibrium solution for the 
gas can be found ( ISutoetalJl9"93) assuming that the baryonic frac- 
tion of the halo is negligible and the gas, with density pg, tempera- 
ture Tg and pressure Pg, follows a polytropic equation of state with 
index Jp, Pg oc pg'' . The equation of hydrostatic equilibrium is: 



dP, _ ^ P9M{< r) 
dr 



(2) 



Here A/(< r) is the halo mass within r (we will call M200 = 
M{< r2oo) the mass within the radius r2oo encompassing an 
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overdensity of 200pcrit), which is given by the NFW profile. 
Under these assumptions the equation can be solved analytically 
dKomatsu & SeljamOOlh : 



Pair) = pgo 



1 - a 1 - 



90 



90 



ln(l + CnfwJ:^) 

Cnfw X 

ln(l + Cni^x) 
Cliffy X 

lri(l + Cnfwa;) 



n i/(7p-i) 



7p/(7p-l) 



; (3) 



Here Pgo, Tgo and Pgo are the gas density, temperature and pressure 
at r = and x — r/r2oo- Furthermore, Cnfw = 7'20o/''s is the 
NFW concentration parameter, while 7^ is the effective polytropic 
index, which determines the shape of the temperature profile. The 
constant a is defined as: 



7p ■ 



1 3 



7p '^O (1 + Cnfw) ln(l + Cnfw) — Cnfw 

The parameter 

3r2oofcBTo 



Vo = 



GnmpM2oo 



(4) 



(5) 



carries information about the central gas temperature To and the 
mass Af2oo, thereby fixing the normalization of the polytropic 
equation of state. In eq. l|5), is the Boltzmann constant, /j, is 
the mean molecular weight (~ 0.58 for a plasma of primordial 
composition), and rUp the proton mass. 

Therefore, pgo is fixed by the constraint on the total gas mass 
Afg within the virial radius. Calling I the integral 



1 



ln(l + t) 
t 



(6) 



(where for simplicity we declare only the dependence on the a ex- 
ponent), we have 



P90 = 



Mg 

47rr§ 



2:(1/(7p - 1)) 



(7) 



In this way, rjo is fixed by the constraint on the total gas thermal 
energy Eg within the virial radius: 



Eg = X 2:(7p/(7p - 1)) ■ 

pmp 



(8) 



The solution of the two conditions must be found numerically by 
iteration. 



2.2 The classical cooling model 

Most SAMs describe the cooling of gas following the model of 
IWhite&Freri3 ( ll99lh . The system is assumed to be in hydrostatic 
equilibrium at the time t = Q. For each mass shell at radius r a 
cooling time can be defined as: 



dlnT 



dt 



3kBTg{r)iimp 
2pg{r)A{Tg{r)) 



(9) 



where A is the cooling function. The left-hand-side of the above 
equation follows from the assumption that dT is computed for 
an isobaric transformation. In the following, we assume, for both 
simulations and analytical models, A(T) to be that tabulated by 



[Sutherland & Dopitj for zero metallicit}^. The cooling ra- 

dius at the time t for the classical model, rc{t), is the defined 
through the relation 



rc{t) ■ tcooi{rc) = t. 



(10) 



In other words, the function rc{t) is the inverse of the function 
tcooi{r). It is then assumed that each shell cools after one cooling 
time. The resulting mass deposition rate reads 



iW'cooi = 4:nr'^pg{rc)^^ . 

at 



(11) 



where a dot denotes a time derivative. 

We emphasize that the classical cooling model imposes that 
a shell of gas cools exactly after one cooling time, computed on 
the initial configuration. The mass cooling rate of equation Jl 11 1, 
predicted by this model, is generally not far from that of the self - 
similar solutions of cooling flows found by iBertschingej ( Il989l) . 
This point will be further discussed in Section 5 and in the Ap- 
pendix. 



2.3 The MORGANA cooling models 

In the classical cooling model, hot gas is assumed to be located 
outside the cooling radius rc- This same assumption is used in the 
MORGANA cooling model; we call this cooling radius tm to distin- 
guish it from the classical one. The equilibrium gas configuration 
is computed as in equations ([3}, but taking into account that no hot 
mass is present within tm', this implies that the integral in equa- 
tion ^ is evaluated from ru / Ts to Cnfw . 

Given the equilibrium profile, the cooling rate of a shell of gas 
of width Ar at a radius r is: 



AAfcooi(r) = 



47rr^Pg(r) Ar 



(12) 



icooi(r) 

where t^ooi is given by equation {9](. The mass deposition rate 
is then computed by integrating the contributions of all the mass 
shells. In performing this integral, we note that the cooling time 
depends on density and temperature, the density dependence be- 
ing always stronger since the temperature profile is much shallower 
than the density profile. The integration in r can then be performed 
by assuming Tg{r) ~ Tg{rM), while taking into full account the 
radial dependence of the density. The resulting total mass deposi- 
tion rate can then be written as 



A^c, 



(13) 



1 - a 1 - 



ln(l + t) 



2/(7p-l) 



r'dt. 



We apply this mass cooling rate starting from t — tcooi(O) (the 
cooling time at r = 0), under the hypothesis that nothing cools 
before that time. The rate of thermal energy loss by cooling is com- 
puted in a similar way: 

3kTg{r cool) 'i-Kvlpgo 



Ecool 



2pmp 

Cntw 



tcool.O 



(14) 



1 - a 1 - 



ln(l + t) 



n 2/{7p-i) 



fdt. 



^ The cooling rate per unit volume should be n^riih; we transform the 
cooling function so that the cooling rate is n^A, where n = rie + rii = 
p I fmip . For sake of simplicity we do not make this explicit in the equation. 
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The lower extreme of integration comes from the assumption that 
the hot and cold phases are always separated by a sharp transition, 
taking place at the cooling radius rm ■ By equating the mass cooled 
in a time interval dt with the mass contained in a shell dr, one 
obtains the evolution of the cooling radius: 



47rpg(rM)rf 



(15) 



The evolution of the system is then followed by numerically in- 
tegrating the equations l ll3t . l ll4t and l ll5t , starting after a time 
icooi(O). The integration is performed with a simple Runge- 
Kutta integrator with adaptive time-step and the gas profile is re- 
computed at each time step. Clearly, as long as the mass and ther- 
mal energy of each cooled shell is removed from the profile, the 
rest of the gas is unperturbed so its profile does not change with 
time. 

The two main assumptions at the base of the MORGANA cool- 
ing model are: (i) all mass shells contribute to cooling according to 
equation l ll2b . (ii) there is anyway a sharp transition from the hot 
to the cold phase, which guarantees the presence of a sharp border 
separating the regions where cooled and hot gas resides. 

Equation l |15t is valid under the assumption that the gas is 
pressure-supported at the cooling radius, which is clearly false in 
general. In particular, as long as the time derivative of the cooling 
radius is much larger than the gas sound speed, the boundary of the 
cooled region propagates so quickly that any gas motion can be ne- 
glected. This holds only at early times and, therefore, for very low 
rM values. On the contrary, the late evolution of cooling is better 
represented by letting the "cooling hole" close at the sound speed. 
We then define a further cooling radius rjvi.ch, whose evolution is 
given by the equation: 



'"M.ch 



47rp3(rM,ch)r2j 1^ 



s(rM,ch) ■ 



(16) 



Here Cs is the sound speed evaluated at the cooling radius. This sim- 
ple change strongly influences the predictions of the cooling model. 
In fact, the evolution of the cooling radius turns out to be remark- 
ably different, with rM.ch < tm; as a consequence, because the gas 
profile is recomputed at each time-step, the mass re-distributes over 
the available volume, at variance with the previous case in which 
the gas profile remains unperturbed beyond rM. We will refer to the 
models of equations l |15t and l |16| l as the "unclosed" and "closed" 
MORGANA cooling models, respectively. 

Equation J16t clearly gives an oversimplified description of 
the process in play, with the implicit assumption that the gas has 
time to settle into a new hydrostatic equilibrium configuration. 
However, due to the negative gas temperature profile, the uncooled 
gas is predicted to become progressively colder, with the result that 
the density and temperature profiles become steeper to keep sat- 
isfying the condition of hydrostatic equilibrium. Eventually, this 
leads to a catastrophic cooling involving all the the gas of the halo. 
A simple and effective way to obtain an acceptable behaviour for 
Mcooi is to suppress the sound speed term of equation l ll6t when 
the specific thermal energy of the hot gas becomes smaller than the 
specific virial energy (—0.5 Uh /Mh, where Uh is the binding en- 
ergy of the DM halo). Since this is admittedly an ad-hoc solution 
to the above problem, our attitude is to consider the "closed" MOR- 
GANA model as an effective model to describe the evolution of the 
cooled mass in DM halos. 

In summary, we have identified three analytic cooling models: 
classical, unclosed MORGANA and closed MORGANA. The main 
difference between the classical and unclosed MORGANA model 



is the following: while the classical model equates the time re- 
quired for a shell to cool to low temperature with the cooling time 
computed on the initial conditions (eq.|9]l, the unclosed MORGANA 
model does not rely on this strong assumption. The main differ- 
ence between the unclosed and closed MORGANA models is that 
the latter attempts a more realistic description of the evolution of 
the cooling radius, taking into account that the cooled gas cannot 
provide pressure support to the inflowing cooling gas. 



3 THE SIMULATIONS 

The simulations that we will discuss in this paper have been ob- 
tained by evolving initial conditions for gas sitting in hydrostatic 
equilibrium within isolated halos whose DM density profile is the 
NEW one. They have been performed using GADGET-2 code, 
a massively parallel Tree-l-SPH code (Springel 20Q5) with fully 
adaptive time-step integration. The version of the code that we 
used adopts an SPH formulation with entropy conserving integra- 
tion and arithmetic symme trization of the hydrodynamical forces 
dSpringel & Hem quist 2002), and includes radiative cooling com- 
puted for a primordial plasma with va nishing metallicity. Th e above 
SPH formulation is also that used bv lYoshida et al.l ( |2002[) in their 
comparison between SAMs and hydrodynamical simulations. It en- 
sures the suppression of spurio us cooling at the interf aces between 
cooled and hot gas (see also iTomatore et al] l2003h . In order to 
follow in detail the trajectories of gas particles in the phase dia- 
gram, while they are undergoing cooling, we have implemented a 
quite conservative criterion of time-stepping, in which the maxi- 
mum time-step allowed for a gas particle is one tenth of its cooling 
time. 

In the following, we will present simulations based on includ- 
ing radiative cooling along with a simple recipe for star formation, 
but excluding any form of energy feedback from supernovae. Only 
in one case, in which we want to study the structure and the evolu- 
tion of the phase diagrams, we turned star formation off. The star 
formation recipe adopted assumes that a collisional gas particle, 
whose equivalent hydrogen number density exceeds uh = 0.1 
cm^"* and with temperature below 3 x 10*K, is instantaneously 
converted into a collisionless "star" particle. The practical advan- 
tage of including star formation is that the simulations becomes 
computationally much faster, since one avoids performing intensive 
SPH computations among cooled high-density gas particles. As we 
shall discuss in the Section 4, we verified that the evolution of the 
mass deposition rate is essentially independent of the introduction 
of star formation. 

Initial conditions for isolated halos have been created by plac- 
ing gas in hydrostatic equilibrium within a DM halo, having the 
NEW density profile (see eq. [T), according to the model given in 
section [27T| To fix the gas therma l energy, we required, as sug- 
gested bv lKomatsu & SeliakI J200lh . that the slopes of the DM and 
gas density profiles be equal at the virial radius. This leads to ther- 
mal energies very sim ilar to 1.2 times the virial energy, as used by 
iMonaco et al.l l l2007h . In order to generate initial conditions, initial 
positions of DM and gas particles are generated by Monte Carlo 
sampling the analytical profiles of eqs. O and (O. To create an 
equilibrium configuration for the DM halo, initial velocities of the 
particl es are assigned according to a local Maxwellian approxima- 
tion ( Hemauistlll993h . where the width of the distribution is given 
by the velocity dispersion of the DM particles, as obtained by solv- 
ing the Jeans' equation. As for the gas particles, their internal en- 
ergy is assigned according to the third of the equations (O. Since 
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Figure 1. Profiles of gas temperature {Top panel), DM density (Bottom left panels) and gas density (Bottom right panel) for tlie non-radiative run of tlie halo 
HI, evaluated at 7 different epochs in the interval [0, 10]trfy„. 



Table 1. Characteristics of the simulated halos. Column 1: halo name; col- 
umn 2: halo mass enclosed within r200\ Column 3: NFW concentration 
parameter; Column 4: redshift used to compute the reference critical den- 
sity; Column 5: value of r200 (in kpc); Columns 6: dynamical time (in Gyr); 
Column 7: central cooling time (in Gyr). 
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the above equations are a hydrostatic solution, gas particles are ini- 
tially assigned zero velocities. 

Each halo has been sampled with 6 x 10 DM particles in- 
side r2oo and an initially equal number of gas particles. The ratio 
between the mass of DM and gas particles is determined by the 
baryon fraction, that we assume to be /bar = 0.19. As for the 
choice of the gravitational softening, it has been chose n to be about 
three t imes larger than the lower limit recommended bv lPower et al.l 
1 I2OO3I) , e ~ 3r2oo /VN200 for a Plummer-equivalent softening, 
where N200 is the number of particles within r2oo ■ The minimum 



value allowed for the SPH smoothing length is assumed to be 0.5 
times the value of the gravitational softening, with SPH computa- 
tions performed using A'ngb = 32 for the number of neighbours. 

To ensure stability of the halos in the absence of cooling, 
density profiles have been sampled with particles out to about 
8r2oo- This also ensures an adequate reservoir of external gas than 
can flow in while cooling removes pressure support in the central 
halo regions. In principle, our controlled numerical experiments 
could have been performed by using a static NFW potential, in- 
stead of sampling the halo with DM particles. However, this pro- 
cedure would not have allowed us to account for any backreac- 
tion of cooling on the DM component. Indeed, it is known from 
cosmological simulations that including gas cooling causes a size- 
able c hange (adiabatic co ntraction) of the structure of the DM halos 
(e.g., iGnedin et al.1l2004l) . 

The characteristics of the simulated halos are summarized in 
Table 1. In this table we also specify the redshift at which the simu- 
lated halo is assumed to stay. Although redshift never explicitly en- 
ters in our simulations of isolated halos, it appears in a indirect way 
when we fix the value of the critical density. Ultimately, increasing 
redshift amounts to take a higher value of pcrit and, therefore, a 
higher halo density for a fixed value of M2oo- In this way, we ex- 
pect that "high-redshift" halos will have a shorter cooling time. In 
the following, we assume the relation between redshift and critical 
density to be that of a cosmological model with Qrn ~ 0.3 and 
Qa = 0.7. 
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Figure 2. Top panels: Density of gas particles for the non-radiative run of tlie lialo HI, as a function of their halo-centric distance, at i = 0, as assigned in the 
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Our reference halo has a mass M200 = 10" M0 (HI in Ta- 
ble 1), typical of a poor galaxy group, with a value of the concen- 
tration parameter Cnfw = 6.3, given by the relation between mass 
and concentration provided by NFW, and r2oo corresponding to the 
value of pcrit al z = Q. Two other halos of the same mass are also 
simulated, which have different values of the concentration param- 
eter, still lying within the scatter in the M2oo-Cnfw relation (H2 and 
H3). We then simulate a smaller (H4) and a larger (H5) halo, with 
M200 = lO^^M© and M200 = IO^^Mq, respectively, so as to 
sample also the scale of elliptical galaxies and of rich clusters. We 
finally consider two halos at z = 1 (H6 and H7) and one halo at 
z = 2 (H8). In all runs we used the same value, 7^ = 1.18, for the 
effective polytropic index. In Table I we also report for each halo 
the values of the dynamical time-scale, which is defined as 



1/2 



^dyn — 



47rGp 



(17) 



and that of the cooling time calculated at the center of the halo as 
in equation 

In order to verify the numerical stability of our results, we 
also performed the following runs for the HI halo: (a) & simula- 
tion with a 10 times larger number of particles and a rescaling of 
the softening, in order to check the effect of mass resolution (HI- 
HR); (b) a. simulation with a 4 times larger number of particles (at 
fixed halo mass) and number of neighbours A'ngb than in the refer- 
ence run, keeping the softening constant for the gravitational force 



(H1-4SPH): because mass resolution is given by A'ngb times par- 
ticle mass, this is kept constant while decreasing the discreteness 
noise in the SPH computation; (c) a simulation with cooling only 
and without star formation (Hl-C); (d) a simulation in which the 
gravitational softening is halved with respect to the reference value 
(Hl-S). 

All the initial conditions have been first evolved for 10 dynam- 
ical times, without cooling. This allowed us to check that tempera- 
ture and density profiles are always stable, thus confirming that the 
initial conditions are indeed quite close to configurations of hydro- 
static equilibrium. An example of this is shown in Figure[T] where 
we plot the profiles of gas density, DM density and temperature for 
the reference halo HI, at different epochs. Although the profiles are 
all remarkably stable, it is worth reminding that there are at least 
two reasons why our initial conditions may not be equilibrium con- 
figurations. Firstly, the gas profiles given by equations ([3) are not 
an exact equilibrium solution because the gas mass is not negligi- 
ble. Secondly, initial particle positions are assigned by performing 
a Monte Carlo sampling of the gas and DM density profiles, while 
the internal energy of gas particles is assigned in a deterministic 
way, by using the third of the equations l|3). The scatter associated 
to the Poissonian sampling of the density profiles, joined with the 
deterministic assignment of the internal energy of SPH particles, 
may well not represent a configuration of stable equilibrium. If this 
is the case, we then expect that the system relaxes to the true min- 
imum energy configuration in a time comparable to its dynamical 
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time-scale. Indeed, Figure|2]shows the scatter plot of individual val- 
ues of density and temperature of all the gas particles, as a function 
of their halo-centric distance, for the initial conditions of Hi (left 
panels) and for a configuration obtained by evolving the system for 
2.5tdyn in the absence of cooling (right panels). The relaxed con- 
figuration shows residual scatter both in density (amounting to 5 
per cent) and temperature (amounting to 15 per cent). The same 
amount of scatter is also found in the high resolution run (Hl-HR), 
while a smaller amount (3 per cent in density, 5 per cent in tem- 
perature) is found in the H1-4SPH simulation, where the density is 
computed averaging over 128 instead of 32 particles. The latter re- 
sult suggests a numerical origin for this scatter, related to the finite 
number of particles used in the SPH computations. 

In order to account for this relaxation, we decided to use as 
initial conditions for the radiative runs the configuration attained by 
evolving the non-radiative runs for 2.5 dynamical times. Once cool- 
ing is turned on, all simulations are let then evolving for 8 central 
cooling times, with the exception of H5, which is run for roughly 
two central cooling times. 



4 RESULTS 

As already emphasized, the main aim of our analysis is to under- 
stand the radiative cooling of the gas in DM halos and to point 
out which one of the cooling models, described in section 2, gives 
results on the evolution of the cooled gas mass which are in best 
agreement with the numerical experiment. To this purpose, most of 
the discussion on how gas cools in simulations will refer to the HI 
halo. Since the evolution of A/cooi is the central result of our anal- 
ysis, we will show it for all the halos and compare the simulation 
results with the predictions of the cooling models. 

In Figure [3] we show the evolution of il/cooi for the reference 
run of the halo HI (i.e. cooling and star formation), and compare 
it with the same run without star formation (Hl-C), with the run 



having 10 times better mass resolution (Hl-HR), four times more 
particles and SPH neighbours (H1-4SPH) and with the run with 
standard mass resolution but with gravitational softening smaller 
by a factor of two (Hl-S). In the runs with cooling and star for- 
mation, A/cool is contributed both by the mass in coUisionless stars 
and by the mass in cold gas particles, which have temperature be- 
low 3 X 10'* K. In the run with cooling only, Mcooi is clearly con- 
tributed only by particles colder than the above temperature limit. 
Figure [3] shows that the evolution of the cooled mass is indepen- 
dent of whether cold and dense particles are treated as coUision- 
les s or SPH particles. T his result, which agrees with that found 
bv lTomatore et al.l j2003l) for cosmological simulations of clusters, 
confirms that using the SPH scheme with explicit entropy con- 
servation and arithme t ic sym metrization of hydrodynamical forces 
ISpringel & Hemguistl ( |2002|) is able to suppress the spurious gas 
cooling which otherwise takes place at the interface between cold 
and hot phases. This result also demonstrates that including star for- 
mation in the radiative runs does not affect our results on Afcooi . As 
for the run with a larger number of neighbours (H1-4SPH), cool- 
ing turns out to start at a later epoch. The reason for this lies in the 
reduced scatter in density in the initial conditions, when a larger 
number of neighbours for the SPH computations is used. In fact, a 
gas particle whose density is scattered upwards with respect to the 
density computed from the profile at its halo-centric distance, has 
a shorter cooling time. As a consequence, the smaller the scatter, 
the lower the probability that a gas particle has cooling time signif- 
icantly shorter than that relative to its radial coordinate. As for the 
high-resolution run (Hl-HR), the larger number of particles pro- 
vides a better sampling of the scattered density distribution. There- 
fore, there is an increasing probability to have a small number of 
gas particles, with exceptionally up-scattered density, which cool 
down at earlier times. Finally, decreasing the gravitational soften- 
ing by a factor 2 (Hl-S run) leads to better resolving the very central 
part of the halo, where gas can more easily cool. This turns into a 
stronger initial transient in the evolution of A/cooi at the onset of 
cooling. Quite reassuringly, despite the differences in the evolution 
of Afcooi between these runs during the first cooling time, they all 
nicely converge after about two cooling times. This demonstrates 
that our numerical description of the evolution of the cooled mass 
is numerically stable after an initial transient. 

In order to probe in detail the behaviour of gas in radiative 
simulations, we have randomly selected five particles in three dis- 
tance intervals from the center (10-20 kpc, 30-40 kpc and 50-60 
kpc) in the initial conditions and we have followed their evolution 
in the Hl-C run (the absence of star formation allows to follow 
the transition from the hot to the cold phase). In Figure |4] we plot 
the evolution of temperature and density as a function of time for 
the selected particles. While flowing toward the halo centre, each 
particle roughly maintains its initial temperature while its density 
progressively increases. Afterwards, in a very short time interval, 
it cools down to ^ lO^K, which corresponds to the temperature 
where the cooling function dies. This means that the transition from 
the hot to the cold phase takes place quite rapidly, thus keeping the 
two phases well separated. 

Figure[5]shows the scatter plot of temperature vs. halo-centric 
distance of the gas particles in the Hl-C run at two output times. 
It is possible to identify a "cooling region" as the spherical shell 
where the drop in temperature takes place. It is interesting to note 
that the size of this region is roughly constant in time and com- 
parable to the softening scale. We have verified that the presence 
of a sharp physical boundary separating hot and cold gas phases 
is robust against numerical resolution, in fact with an even sharper 



8 Viola et al. 



le+08 
le+07 
^ le+06 
H 100000 
10000 



i 1 1 1 1 1 — 1 — 1 — 1 1 i 


1 1 1 1 — 1 — 1 — 1 — 1 1 : 






ri 










! 1 

r ■ 1 

- ; 1 




' i 
/ / i 
' / J 




1 ' 






1 





le-21 
le-22 
le-23 
le-24 
\ le-25 



a. 



le+08 
le+07 
^ le+06 
H 100000 
10000 



E 1 1 1 1 1 1 — 1 — 1 1 : 


1 1 1 1 — 1 — 1 — 1 — 1 1 : 








i :l ; 
i ii ; 


r i i 
■ i i 






■' v/ ;' - 








/ .•■// / 


1 1 1 1 1 1 1 1 1 





le-21 
le-22 
le-23 
le-24 
le-25 



a. 



le+08 
le+07 
g le+06 
H 100000 
10000 




le+09 
t(yr) 



le+09 
t(yr) 





le-21 






le-22 


cn 




le-23 


B 




le-24 






le-25 


a. 



Figure 4. Evolution of temperature (left panels) and density (right panels j /or three sets of five particles each, selected within three different radial shells in the 
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boundary at higher resolution, but its size does depend on resolu- 
tion. In the H 1 -HR run, who has roughly a two times smaller soften- 
ing, the size of the cooling region is more than 50 per cent smaller. 
Therefore, while our simulations provide a numerically convergent 
result on the evolution of A/cooi, they do not provide a similarly 
convergent result on the size of the cooling region. 

In Figure|6]we show the density profile of DM and the density, 
pressure and temperature profiles of the hot gas in the simulation 
at several times, normalized to the corresponding profiles evaluated 
for the initial conditions. Pressure increases in the inner part of the 
halo for a few dynamical times, to saturate later to values that peak 



at a factor of six times higher than in the initial conditions, just out- 
side the cooling region. This increase in pressure is mainly driven 
by a comparable increase in gas density, while the temperature pro- 
file is much more stable or even slightly decreasing near the cooling 
region. The density increase on the other hand is partly caused by 
the adiabatic contraction of the DM halo, but because the increase 
in DM density, though significant, is smaller than that in gas den- 
sity, then adiabatic contraction gives only a minor contribution. As 
we shall discuss in the following, the increase in gas density en- 
hances radiative losses as the gas approaches the cooling region. 
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and this turns into a shortening the cooHng time of the inflowing 
gas particles, with respect to the predictions of the classical model. 

In summary, the following conclusions can be drawn from our 
analysis of gas cooling in simulations: 

(i) the drop in temperature when gas particles pass from the hot to 
the cold phase is quite rapid; 

(ii) this drop in temperature takes place in a spherical shell with a 
rather sharp boundary, the cooling region, which separates the inner 
and outer regions dominated by cold and hot particles, respectively; 

(iii) density and pressure increase with time just beyond the cooling 
region, and this increase is not driven by adiabatic contraction of 
the halo. 

In the light of these results, the question then arises as to 
whether the cooling recipes, described in section 2, able to repro- 
duce the rate of mass cooling found in the simulations. 

In order to answer this question, we first address the issue con- 
cerning the tuning of the initial conditions. As mentioned in section 
3, the gas in the simulations relaxes to a minimum energy configu- 
ration, so that the parameters of the gas profile after 2.5 dynamical 
times, when the cooling is turned on, may in principle differ from 
the ones used to generate the initial conditions. This is a crucial 
point because in order to make a reliable comparison between an- 
alytic cooling models and numerical simulations, we must be sure 
that the initial conditions are the same for both. Owing to the sta- 
bility of the profiles shown in Figure [T] we expect this effect to be 
small. As we shall see, even small differences in the initial pro- 
files leads to an appreciable change in the resulting evolution of 
Mcooi • In the model that we used to generate the initial hydrostatic 
configurations (equations [3}, the parameters that determine the gas 
density and the temperature profiles are the halo gas mass Mg, the 
effective polytropic index 7p, and the central temperature of the 
gas (in units of the virial one) r]. While holding the mass fixed, 
we have considered 900 pairs of values for the parameters (7p, rf), 
varying both the polytropic index and the energy factor in the range 
[1.1 : 1.4]. We have calculated for each pair the theoretical density 



and temperature profiles according to equations JS}, and compared 
them with the profiles from the initial conditions. In particular, we 
calculated for each radius the root mean square difference between 
the density and temperature profiles of the hydrostatic model and 
those of the initial conditions, and imposed the maximum differ- 
ence to be smaller than 10 per cent. Then, for each cooling model 
and unless otherwise stated, we show predictions relative to all the 
profiles that were selected by the procedure described above. This 
allows us to keep control on any uncertainty of the initial profiles 
which are used as input to the cooling models. 

As a main term of comparison between analytic models and 
simulations, we use the evolution of the cooled mass fraction. In 
Figure |7] we plot the evolution of Mcooi for all the eight simulated 
halos from simulations and the predictions of the different cooling 
models. As for the latter, each shaded area represent the envelope 
of each model prediction for all the profiles which provide a good 
fit to the initial conditions. 

From these plots we infer the following conclusions. 

(i) No gas is allowed to cool down to low temperatures before one 
central cooling time has elapsed, both in the classical and in the 
MORGANA models. Afterwards, cooling takes place abruptly, giv- 
ing rise to a sort of "burst" of star formation. This behaviour is 
consistent with simulation results, at least when the scatter in den- 
sity and temperature, which are present in the initial conditions, is 
reduced. 

(ii) The classical cooling model (red area) always underpredicts the 
value of Afcooi- This underestimate is more severe at earlier times, 
and remains quite substantial, a factor of 2-3, even in the most 
evolved configurations. This effect is generally stronger for the ha- 
los having a lower concentration and/or higher mass, i.e. having 
longer cooling times (see Table 1). 

(iii) The unclosed MORGANA model (magenta area) follows in a 
much closer way the cooled mass at early times and the fit is very 
good in most cases at epochs between one and a few central cooling 
times. In the most evolved configurations, the model underestimate 
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Figure 6. Profiles of DM density (upper left panel), gas density (upper light), pressure (lower left) and temperature (lower right) at different epochs for the HI 
simulation with only cooling (Hl-C). All profiles are normalized to their initial value. 



the cooled mass, although the underestimate is sensibly reduced 
with respect to the classical model. 

(iv) The closed MORGANA model behaves very similarly to the un- 
closed one for a few central cooling times. At later times, it predicts 
a larger value of i\fcooi, providing a generally good fit to the simu- 
lation results for the most evolved configurations. 



5 DISCUSSION 

The better performance of the unclosed MORGANA model, with re- 
spect to the classical model, in reproducing the evolution of the 
cooled gas fraction from the simulations can be ascribed to the fol- 
lowing facts. 

As explained in Section 2.4, the classical model relies on the 
strong hypothesis that each gas shell cools in exactly one cooling 



time, with this cooling time computed on the initial conditions. 
While this hypothesis is valid when the first gas particles cool, it 
breaks down soon later. This behaviour is indeed not unexpected. 
The evolution of a mass element in presence of cooling and adia- 
batic compression is given by the following equation: 

f^T[~-^ + ^n, (18) 

\ Icool op/ 

where p{t) and T{t) describe the evolution of density and temper- 
ature of the gas element, and tcooi is the cooling time computed on 
the actual density and temperature (not on the initial conditions). 
Under the assumptions that the temperature dependence of the 
cooling time can be neglected, so that tcooi(t) ~ tco{p{t)/ po)~^ , 
and that pressure is constant during the evolution, it is easy to solve 
equation l |18t and find that the time ttot required for the mass el- 
ement to cool to T = coincides with tco- Therefore, the basic 
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assumption of the classical model is satisfied in the case in which 
gas particles cool at constant pressure. 

On the other hand, Figure|4] shows that gas particles take most 
time to flow from their initial position to the cooling region. Since 
cooling takes place in a pressure-supported way during this time, 
their radiative losses are balanced by adiabatic compression. As 



a result, the temperature of the gas particles has a slow evolution 
while density increases significantly as they move towards the cool- 
ing region. This results in shallow and stable temperature profiles, 
while density profiles become progressively steeper (see Figure|6]l. 
The density increase turns into enhanced radiative losses, thereby 
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making the total cooling time, ttot, shorter than the cooling time, 
tco, computed on the initial conditions. 

The main assumption of the classical model is then clearly in- 
validated by our simulations. Going ba ck to the original proposal 
of this model bv lWhite & FrenkI h99lh . the main justification was 
that the model i s roughly consistent with the exact self-similar so- 
lutions found bv lBertschin ger (1989). However, Bertschinger's self 
similar solutions give cooling flows that are equal to: 



Mcooi = 47rr-^pg(rc)^^ X , 
at 



(19) 



where the constant /lo depends on the initial profile and can take 
values ranging from ~ 0.1 to ~ 2.5 (see tables 1 and 2 in 
Bertchinger's paper). The classical model is recovered in the case 
Ho = 1. In the Appendix we will consider the simple case of an 
isothermal halo with power-law density profile pg oc r~^. In this 
case the unclosed MORGANA model gives ttot ~ O.SicO and a 
mass-deposition rate, A/cooi, which is equal to \/2 times that of the 
classical model. On the other hand, Bertschinger's solutions give 
cooling flows higher than the classical value by a factor ranging 
from 1.304 to 1.190, depending on the shape of the cooling func- 
tion, thus implying total cooling times ttot shorter than tco by a fac- 
tors ranging from 0.59 to 0.70. Both MORGANA and Bertschinger's 
self-similar solutions predict that shallower power-law profiles give 
stronger cooling flows and shorter total cooling times. Therefore, 
while the unclosed MORGANA model agrees with the exact self- 
similar solution to within a few tens per cent, the assumption that 
ttot ~ tco is generally not correct. 

More in detail, the difference between classical and unclosed 
MORGANA models at the onset of cooling is the consequence of 
the flatness of the density profile in the inner regions of the gas. 
This can be shown by finding exact solutions of the classical 
and MORGANA models in the case of power-law density profiles, 
Pg(r) oc r~", for an isothermal gas. These calculations are re- 
ported in the Appendix and the results can be summarized as fol- 
lows, (i) Both the unclosed MORGANA model and the classical one 
predict a self-similar cooling flow for a > 3/2, thus implying that 
the corresponding mass deposition rates are proportional to each 
other, (ii) The classical and unclosed MORGANA models agree if 
a — 3; shallower profiles lead the unclosed MORGANA model to 
predict shorter total cooling times and higher cooling rates, (iii) For 
profiles flatter than q = 3/2 the mass cooling flow of the unclosed 
MORGANA model is dominated by the external regions. The solu- 
tion is not proportional to the classical one, the cooling flow is not 
self-similar and is roughly constant. Clearly, such a shallow profile 
will hold in the inner region of a realistic halo. 

As a conclusion, the strict validity of the classical model is 
limited to specific profiles and to self-similar flows. On the other 
hand, the unclosed MORGANA model, which relaxes the assump- 
tion on the total cooling time of a gas shell, better reproduces the 
stronger flows found in our simulations of isolated halos, which 
takes place when the Lagrangian cooling radius sweeps the flat part 
of the density profile. 

The closed MORGANA model further improves the agreement 
with the simulations by increasing the fraction of cooled mass. The 
main reason for this increase is that, being always tm.cIi < '"m, the 
smaller value of the cooling radius leads to an increase of the den- 
sity of the cooling shell, simply because the hot gas is allowed to 
stay at smaller radii. This implies still shorter cooling times and en- 
hanced cooling flows. Another prediction of the closed MORGANA 
model is that the cooling radius rM.ch is stable after a quick tran- 
sient. This prediction is in qualitative agreement with the results of 



the simulations (see Figure [5j the dotted lines give the position of 
'"M.ch). However, the size of the cooling region in the simulation 
is affected by resolution, so this comparison cannot be pushed to 
a quantitative level. The validity of this model breaks as soon as 
the energy of the uncooled gas drops below the virial value. In this 
case we still obtain a good match of the cooled mass fraction by 
simply dropping the sound speed term in equation il6\ . which is 
responsible for the shrinking of the cooling hole. For this reason, 
we consider the closed MORGANA model as an effective model, in 
that it takes into account the shrinking of the cooling region caused 
by the pressure from the hot gas just outside this region, without 
however providing a formally rigorous description for this effect. 

Using the predicted rough constancy of the cooling flow when 
the central, shallow part of the gas profile is cooling, in Appendix 
B we show that it is possible to give a very simple and remarkably 
accurate prediction of the cooled mass as the result of a constant 
flow which is given by a simple analytic formula, valid up to ~ 10 
central cooling times. 



6 CONCLUSIONS 

We have presented a detailed analysis of cooling of hot gas in DM 
halos, comparing the predictions of semi-analytic models with the 
results of controlled numerical experiments of isolated NFW halos 
with hot gas in hydrostatic equilibrium. Simulations have been per- 
formed spanning a range of masses (from galaxy- to cluster-sized 
halos), concentrations and redshift (from to 2). Smaller halos at 
higher redshift have not been simulated because the validity of the 
assumption of a hydrostatic atmosphere is doubtful when the cool- 
ing time is much shorter than the dynamical time. 

We have consi dered the "classical" cooling model of 
IWhite&Fren3(ll99ll), used in mo st SAMs, and the model recently 
proposed by Monaco et alj llOOf) within the MORGANA code for 
the evolution of galaxies and AGNs. The main features of these 
models can be summarized as follows. The density and pressure 
profiles of the gas are computed by solving the equ ation of hydro- 
static equilibrium in an NFW halo ( ISutoet alJll998h . The classical 
cooling model assumes that each mass shell cools to low tempera- 
ture exactly after one cooling time tcooi{r), computed on the initial 
conditions. The cooling radius rc is then the inverse of the tcooi{r) 
function, and the cooled fraction is the fraction of gas mass within 
rc. The "unclosed MORGANA" cooling model computes the cool- 
ing rate of each mass shell, then integrates over the contribution of 
all mass shells and follows the evolution of the cooling radius as- 
suming that the transition from hot to cold phases is quick enough 
so that a sharp border in the density profile of hot gas is always 
present. This determines the evolution of the cooling radius rm- 
Moreover, to mimic the closure of the "cooling hole" due to the lack 
of pressure support at tm, the cooling radius (now called tm.cIi) is 
induced to close at the sound speed. This defines the "closed MOR- 
GANA" model. 

Our main results can be summarized as follows, 
(i) The classical cooling model systematically underestimates the 
fractions of cooled mass. After about two central cooling times, 
they are predicted to be about one order of magnitude smaller than 
those found in simulations. Although this difference decreases with 
time, after 8 central cooling times, when simulations are stopped, 
the difference still amounts to a factor 2-3. This disagreement is as- 
cribed to the lack of validity of the assumption that each mass shell 
takes one cooling time, computed on the initial conditions, to cool 
to low temperature. Seen from the point of view of a mass element. 
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the time required by it to cool to low temperature is shorter than 
the initial cooling time when density increases and temperature is 
constant during cooling. This is what happens to gas particles in 
the simulations: they take most of time to travel from their initial 
position toward the cooling region, at roughly constant tempera- 
ture and increasing density. The disagreement is stronger when the 
cooling gas comes from the shallow central region, in which case 
the cooling flow is markedly not self-similar. 

(ii) The unclosed MORGANA model gives a much better fit of the 
cooled mass fraction. This is mostly due to the relaxation of the 
assumption on the cooling time, mentioned in point (i). This model 
correctly predicts cooling flows which are stronger than the clas- 
sical model, by a larger amount for flatter gas density profiles. In 
the Appendix we show that the solution is not self-similar if the 
slope of the density profile is shallower than r~'^^^. In this case 
cooling is not dominated by the the shells just beyond the cooling 
radius but the whole region for which the density profile is shallow 
contributes. 

(iii) The closed MORGANA model further improves the fit to the 
simulation results on the evolution of the cooled mass fraction, giv- 
ing accurate results to within 20-50 per cent in all the considered 
cases, after about 8 central cooling times. This agreement is a good 
reward for the increase of physical motivation of this model, ob- 
tained at the modest price of letting the cooling radius close at the 
local sound speed. However, the closure of the cooling radius must 
be halted at later times for the model to give realistic results. In 
general, we consider the closed MORGANA as a successful effec- 
tive model of cooling, rather than as a rigorous physical model. 

(iv) The cooling flow is well approximated by a constant flow, for 
which we give a fitting formula in Appendix B, and which is valid 
up to ~ 10 central cooling times. 

In the context of models of galaxy formation, cooling of hot 
virialized gas is the starting point for all the astrophysical processes 
involved in the formation of stars (and supermassive black holes) 
and their feedback on the interstellar and intra-cluster media. We 
find that the classical model, used in most SAMs, leads to a sig- 
nificant underestimate of the cooled mass at early times. This re- 
sult is apparently at variance with previous claims, discussed in 
the Introduction, of an agr eement of models and simulations in 
predicting the cooled mass, jBenson et al.l l200ll : iHellv et alj|2003l : 
ICattaneo et"Zll2007l : I Yoshidaet al. 200Z) . Given the much higher 
complexity of the cosmological initial conditions used in such anal- 
yses, it is rather difficult to perform a direct comparison with the 
results of our simulations of isolated halos. Here we only want to 
stress the advantage of performing simple and controlled numerical 
tests in order to study in detail how the process of gas cooling takes 
place. 

It is well possible that the different behaviour of simulations 
and the classical cooling model is less apparent when the more 
complex cosmological evolution is considered. However, there is 
no doubt that the results of our analysis are quite relevant for the 
compari son between SAM p redictions and observations. For in- 
stance, jpontanot et al.l 120071) have recently shown that the MOR- 
GANA model of galaxy formation is able to reproduce the observed 
number counts of sources in the sub-mm band by using the standard 
Initial Mass Function (IMF) by Salpeter ( 1955), without any need 
to resort to a top-heavier IMF CBaugh,2006,) . As argued by these 
authors, the bulk of starbursts are driven by massive cooling flows, 
so this difference is mostly due to the different cooling models. 



ACKNOWLEDGMENTS 

We wish to thank Volker Springel for having provided us with the 
non-public version of GADGET-2 and Fabio Fontanot, Ian Mc- 
Carthy, Richard Bower and Klaus Dolag for many enlightening 
discussions. The simulations have been realized using the super- 
computing facilities at the "Centro Interuniversitario del Nord-Est 
per il Calcolo Elettronico" (CINECA, Bologna), with CPU time as- 
signed thanks to an INAF-CINECA grant and to an agreement be- 
tween CINECA and the University of Trieste. This work has been 
partially supported by the INFN PD-5 1 grant and by a ASI/INAF 
grant for the support of theory activity. 



References 

Baugh C. M., 2006, Reports of Progress in Physics, 69, 3101 
Benson A. J., Pearce F. R., Frenk C. S., Baugh C. M., Jenkins A., 

2001, MNRAS, 320, 261 
Bertschinger E., 1989, ApJ, 340, 666 

Cattaneo A., Blaizot J., Weinberg D. H., Keres D., Colombi S., 
Dave R., Devriendt J., Guiderdoni B., Katz N., 2007, MNRAS, 
377, 63 

Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 
319, 168 

Fontanot F, Monaco P., Silva L., Grazian A., 2007, MNRAS, ac- 
cepted, 0, 

Gnedin O. Y., Ki-avtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 
616, 16 

Helly J. C, Cole S., Frenk C. S., Baugh C. M., Benson A., Lacey 

C, Pearce F R., 2003, MNRAS, 338, 913 
Hernquist L., 1993, ApJS, 86, 389 

Kauffmann G., White S. D. M., Guiderdoni B., 1993, MNRAS, 
264, 201 

Komatsu E., Seljak U., 2001, MNRAS, 327, 1353 
Menci N., Fontana A., Giallongo E., Salimbeni S., 2005, ApJ, 
632, 49 

Monaco P, Fontanot F, Taffoni G., 2007, MNRAS, 375, 1 189 

Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557 

Nagamine K., Cen R., Hernquist L., Ostriker J. P., Springel V., 

2004, ApJ, 610, 45 
Navarro J. E, Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 
Power C, Navarro J. E, Jenkins A., Frenk C. S., White S. D. M., 

Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14 
Romeo A. D., Portinari L., Sommer-Larsen J., 2005, MNRAS, 

361, 983 

Salpeter E. E., 1955, ApJ, 121, 161 

Saro A., Borgani S., Tornatore L., Dolag K., Murante G., Biviano 

A., Calura F, Chariot S., 2006, MNRAS, 373, 397 
Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087 
Springel V., 2005, MNRAS, 364, 1105 

Springel V., Frenk C. S., White S. D. M., 2006, Nature, 440, 1 137 
Springel V., Hernquist L., 2002, MNRAS, 333, 649 
Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253 
Suto Y, Sasaki S., Makino N., 1998, ApJ, 509, 544 
Tornatore L., Borgani S., Springel V., Matteucci F, Menci N., Mu- 
rante G., 2003, MNRAS, 342, 1025 
White S. D. M., Frenk C. S., 1991, ApJ, 379, 52 
Yoshida N., Stoehr E, Springel V., White S. D. M., 2002, MN- 
RAS, 335, 762 



14 Viola et al. 



APPENDIX A: MODEL SOLUTIONS FOR POWER-LAW 
PROFILES 

In this Appendix we will discuss the cases in which the un- 
closed MORGANA cooling model provides self-similar solutions 
and how these solutions compare to the exact self-similar solu- 
tio ns by ISertsc hinger (1989) and to the classical cooling model 
by IWhite & Frenk a991.) . To this purpose, we analytically solve 
equation ( I13t for the evolution of the cooled gas mass Mcooi and 
equation l |15l l for the evolution of the cooling radius tm in the case 
of an isothermal power-law density profile: 



Pair) = P9V 



( — ) 

V r2oo / 



(Al) 



Here pgv is the density at the virial radius r2oo, while temperature 
is fixed to the virial one. In the following we will make the approx- 
imation that we can neglect the dependence of the cooling time on 
temperature: 



(A2) 



/ ^ 

tcV I 
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where tcv is the initial cooling time of the gas at the virial radius. 
Let Tm be the cooling radius of the unclosed Morgana model. The 
mass cooling flow (eauationll3t is: 

-2a + 2 



tcV 
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dx , 



where x = r/r2oo. If a ^ 3/2 the solution is: 
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If Q > 3/2 and r-M ^ ''200 then the equation for the Morgana 
cooling radius vm becomes (equation [Tsj: 



^200 



whose solution is 
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Substituting this solution into equation l lA4b and neglecting the — 1 
term in parentheses we obtain 



47rr: 



20oPaV 



tc 



2a - 3 



\2a-3tcv) 



-(2a-i)/o 



(A7) 



An analogous computation can be performed for the classical 
cooling model. From the definition of equation l llOt for the classical 
cooling radius, we obtain: 



rc{t) = '■200 (7^) 

\tc.V / 

which gives 



l/a 



M,, 



atcv 



(3-2a)/a 



(A8) 



(A9) 



for the mass deposition rate. Therefore, the comparison of the two 
models leads to 



V2a-3y 



{3-a)/a 
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In the case of a = 2 we obtain Mcooi = \/2M:iassic- Since the 
mass deposition rate from the MORGANA model is proportional to 
the classical one, cooling is self-similar and can be compared to the 



solution o f lBertschingeil ( ll989l) . The latter depends on the (power- 
law) temperature dependence of the cooling function, A cx . To 
have a cooling time independent of temperature, as assumed in the 
MORGANA model, we should compare to the solution for A = 1, 
which is not given by Bertschinger because it is of no physical rel- 
evance. Going from A = —1 to A = 0.5 (the values considered 
by Bertschinger) the ratio between the values of Mcooi predicted 
by the self-similar solution and by the classical model grows from 
1.190 to 1.304. This compares well to our prediction of \/2- We 
then conclude that the unclosed MORGANA model is roughly con- 
sistent with Bertschinger's exact self-similar solution within a few 
tens per cent in the case of an isothermal gas density profile with 
a = 2. 

The total cooling time ttot, i.e. the time required by a shell 
to cool to very small temperature, is defined as the inverse of the 
ruit) function. The unclosed MORGANA model gives: 



2a 



-tcV 



2a 



-tcooi(r) . 
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V ^200 / 

while the classical model assumes ttot = tcooi{r). Then, complete 
cooling is predicted by the MORGANA model to take less than one 
cooling time (half of it if a = 2) whenever a < 3. Clearly, the clas- 
sical and unclosed MORGANA model give the same cooling flow 
and total cooling time for an isothermal profile with a = 3. 

If a < 3/2, the —1 term in equation iA4l becomes dominant. 
As a consequence, the cooling flow is not proportional to the classi- 
cal one, so the solution is not self-similar, and cannot be compared 
to Bertschinger's solutions. In fact, in a shallow profiles cooling is 
dominated by the external regions and the flow is roughly constant 
as a first approximation. However, in realistic cases a shallow pro- 
files will be valid only within some reference scale radius rrcf , and 
cooling will proceed very quickly until the cooling radius tm has 
reached r^. 



APPENDIX B: A SIMPLE ANALYTIC EXPRESSION FOR 
THE COOLING FLOW 

The last paragraph of Appendix A shows that, as long as the cool- 
ing radius sweeps the region where the gas density profile is shal- 
lower than r^^''^, the cooling flow is roughly constant. Indeed, 
we find that both simulations and the MORGANA results can be 
fit with a constant cooling flow. To find an analytic approximation 
for it, we start from computing the reference radius r^of at which 
dlnpj/dlnr = —3/2. The function Inpg (equation|3j is Taylor- 
expanded around : 



Pgirs) 
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The reference radius is then: 

—3/2 — dlnpg/dlnr(rs) 



Trcf ^ Ts exp 



d^ In pg/dlnr2(rs) 
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Figure Bl. Simple approximation to tlie cooled mass fraction obtained by integrating in time the constant cooling flow of equation IB 1 1 1 shown for the HI 
and H6 simulations (the others show a similar degree of agreement). Dotted (magenta), dashed (blue) and dot-dashed (green) lines give the predictions of the 
classical, unclosed and closed MORGANA models. The thick continuous line gives the simple analytic fit. For sake of completeness, the models are obtained 
using 7p = 1.20 and thermal energy equal to 1.18 times the virial energy (HI), and 7p = 1.21 and energy equal to 1.15 the virial energy (H6). These values 
are at the centre of the intervals used for the models shown in figure|2] 



A guess for the reference cooling flow can then be obtained as 
M ~ 47rri?(,fPg(r-i.of )/tcooi(''rof )- This guess has been compared to 
the value of the cooling flow in the MORGANA unclosed and closed 
models, averaged over the time interval from 1 to 3 central cool- 
ing times and over the two models. We find systematic differences 
that are removed by adopting a correction function /{Mh, Cnfw) 
of halo mass and concentration. Let's call: 



cooled mass does not overshoot the gas mass. However, a precise 
modeling of these late phases is immaterial, as feedback from star 
formation and AGNs would clearly regulate the dynamics of cool- 
ing, so any reasonable truncation would work equally well. 



A 



1.95 



1.41 X 10' 



B = max <i 0.09 
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Mh 
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Me 



,0.06 



1.26 X 10'^ 
0.043 {B - 0.085) + 0.0014 
then the correction function / is found to be: 

f{MH,Cnl^) =A+ (Cnf„ - 10)5 

if Cnfw <= 10, and 

/(Mh, Cnfw) =A+ (Cnfw - 10)5 + (Cnfw - 10)^C 

if Cnfw > 10. The approximated cooling flow is then: 
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The dependence on cosmology is hidden in the dependence on red- 
shift and concentration of rrcf and the function /{Mh, Cnfw)- 

This simple analytical model gives a very good fit to all the 
simulations shown in this paper, at least at early times. For sake 
of brevity, we show here (figure IbH only the reference HI case 
and the worst-case H6 simulation. In the latter case the fraction of 
cooled mass approaches unity, and the approximation of a constant 
cooling flow breaks after ^ 10 central cooling times. It is very 
simple to truncate the constant cooling flow by imposing that the 



